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While linear mixed modeling methods are foundational concepts 
introduced in any statistical education, adequate general methods 
for interval estimation involving models with more than a few vari- 
ance components are lacking, especially in the unbalanced setting. 
Generalized fiducial inference provides a possible framework that ac- 
commodates this absence of methodology. Under the fabric of gener- 
alized fiducial inference along with sequential Monte Carlo methods, 
we present an approach for interval estimation for both balanced 
and unbalanced Gaussian linear mixed models. We compare the pro- 
posed method to classical and Bayesian results in the literature in 
a simulation study of two-fold nested models and two-factor crossed 
designs with an interaction term. The proposed method is found to 
be competitive or better when evaluated based on frequentist crite- 
ria of empirical coverage and average length of confidence intervals 
for small sample sizes. A MATLAB implementation of the proposed 
algorithm is available from the authors. 

1. Introduction. Inference on parameters of normal linear mixed models 
has an extensive history; see Khuri and Sahai (1985) for a survey of vari- 
ance component methodology, or Chapter 2 of Searle, Casella and McCul- 
loch (1992) for a summary. There are many inference methods for variance 
components such as ANOVA-based methods [Burdick and Graybill (1992), 
Hernandez, Burdick and Birch (1992), Hernandez and Burdick (1993), Je- 
yaratnam and Graybill (1980)], maximum likelihood estimation (MLE) and 
restricted maximum likelihood (REML) [Hartley and Rao (1967), Searle, 
Casella and McCulloch (1992)] along with Bayesian methods [Gelman (2006), 
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Gelman et al. (2004), Wolfinger and Kass (2000)]. Many of the ANOVA- 
based methods become quite complex with complicated models (e.g., due 
to nesting or crossing data structures), and are not guaranteed to perform 
adequately when the designs become unbalanced. When the model design 
is not balanced, the decomposition of the sum-of-squares for ANOVA-based 
methods is not generally unique, chi-squared or independent. Furthermore, 
"exact" ANOVA-based confidence intervals are typically for linear combi- 
nations of variance components, but not the individual variance compo- 
nents even for simple models [Burdick and Graybill (1992), pages 67 and 68, 
Jiang (2007)]; however, approximate intervals do exist. With notable opti- 
mality properties for point estimation, MLE and REML methods are less 
useful when it comes to confidence interval estimation for small samples be- 
cause the asymptotic REML-based confidence intervals tend to have lower 
than stated empirical coverage [Burch (2011), Burdick and Graybill (1992), 
Searle, Casella and McCulloch (1992)]. Bayesian methods, in particular, hi- 
erarchical modeling, can be an effective resolution to complicated models, 
but the delicate question of selecting appropriate prior distributions must 
be addressed. 

There are numerous applications of normal linear mixed models related 
to topics such as animal breeding studies [Burch and Iyer (1997), E, Hannig 
and Iyer (2008)], multilevel studies [O'Connell and McCoach (2008)] and 
longitudinal studies [Laird and Ware (1982)]. Many implemented methods 
do not go beyond two variance components or are designed for a very spe- 
cific setting. We propose a solution based on generalized fiducial inference 
that easily allows for inference beyond two variance components and for the 
general normal linear mixed model settings. 

The proposed generalized fiducial approach is designed specifically for in- 
terval data (e.g., due to the measuring instrument's resolution, rounding for 
storage on a computer or bid-ask spread in financial data). There are several 
reasons to consider interval data. First, there are many examples where it 
is critical [or required per the regulations outlined in GUM (1995)] to incor- 
porate all known sources of uncertainty [Elster (2000), Frenkel and Kirkup 

(2005) , Hannig, Iyer and Wang (2007), Lira and Woger (1997), Taraldsen 

(2006) , Willink (2007)]. Second, our simulation results for the normal lin- 
ear mixed model show that even when considering this extra source of un- 
certainty, the proposed method is competitive or better than classical and 
Bayesian methods that assume the data are exact; see Section 3. [The pro- 
posed method is also appropriate for noninterval, or standard, data simply 
by artificially discretizing the observation space into a fixed grid of nar- 
row intervals; Hannig (2012) proves that as the interval width decreases to 
zero, the generalized fiducial distribution converges to the generalized fidu- 
cial distribution for exact data.] Finally, on a purely philosophical level, all 
continuous data has some degree of uncertainty, as noted previously, due to 
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the resolution of the measuring instrument or truncation for storage on a 
computer. Note that we are not suggesting that all methods should incorpo- 
rate this known uncertainty; however, we were able to appeal to this known 
uncertainty for the computational aspect of the proposed method. 
A general form of a normal linear mixed model is 

(1.1) Y = Xp + VZ + e, 

where Y is an n x 1 vector of data, X is a known n x p fixed effects design 
matrix, (3 is a p x 1 vector of unknown fixed effects, VZ = Yll=x ViZi, where 
Zi is a vector of effects representing each level of random effect i such that 
E{Zi) = and var(Zj) = G, Vj is the known design matrix for random effect 
i, and e is an n x 1 vector representing the error and E(e) = and var(e) = R 
[Jiang (2007)]. Note that there are r total random components in this model, 
and covariance matrices G and R contain unknown parameters known as 
variance components. It is often assumed that Z and e are independent and 
normally distributed. Additional assumptions on this model for the proposed 
method are addressed in Section 2. 

Except where noted, the notational convention that will be used for ma- 
trices, vectors, and single values will be, respectively, bold and capital letters 
for matrices, capital letters for vectors and lowercase letters for single values 
(e.g., A, A, a). 

The focus of this paper is the construction of confidence intervals for the 
unknown parameters of (1.1), with emphasis on the variance components 
of matrices G and R. In this paper, inferences are derived from the gen- 
eralized fiducial distributions of the unknown parameters, and we propose 
a sequential Monte Carlo (SMC) algorithm to obtain these samples. Like a 
Bayesian posterior, this procedure produces a distribution on the parame- 
ter space, but does so without assuming a prior distribution. We evaluate 
the quality of the simulated generalized fiducial distribution based on the 
quality of the confidence intervals — a concept analogous to confidence distri- 
butions [Schweder and Hjort (2002), Xie, Singh and Strawderman (2011)]. 
We begin by introducing the two main techniques of the proposed method: 
generalized fiducial inference and SMC methods. Then we introduce the pro- 
posed method, and state and prove a theorem concluding the convergence 
of the algorithm. To demonstrate small sample performance, we perform 
a simulation study on two different types of models (unbalanced two-fold 
nested models and two-factor crossed with interaction models), and include 
a real-data application for the two-fold nested model. We finish with con- 
cluding remarks. Additional information can be found in the supplemental 
document [Cisewski and Hannig (2012)]. 

1.1. Generalized fiducial inference. Fiducial inference was introduced by 
Fisher (1930) to rectify what he saw as a weakness in the Bayesian philoso- 
phy, where a prior distribution is assumed without sufficient prior knowledge. 
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While Fisher made several attempts at justifying his definition of fiducial 
inference [Fisher (1933, 1935)], it was not fully developed. Fiducial infer- 
ence fell into disrepute when it was discovered that some of the properties 
Fisher claimed did not hold [Lindley (1958), Zabell (1992)]. Efforts were 
made to revitalize fiducial inference by drawing connections to other areas 
such as Fraser's structural inference [Fraser (1961a, 1961b, 1966, 1968)], and 
more recently Hannig, Iyer and Patterson (2006) connect Fisher's fiducial 
inference to generalized inference introduced in Tsui and Weerahandi (1989) 
and Weerahandi (1993). In this paper, we propose a method for inference 
on parameters of normal linear mixed models using the ideas of generalized 
fiducial inference. 

The main idea of fiducial inference is a transference of randomness from 
the model space to the parameter space. A thorough introduction to gener- 
alized fiducial inference can be found in Hannig (2009), but here we consider 
a simple example. Let y be a realization of a random variable Y ~ N(/j,, 1) 
[where N(/j,,l) represents a normal distribution with unknown mean \x and 
standard deviation 1]. The random variable Y can be represented as Y = 
H + Z where Z ~ N(0, 1). Given the observed value y, the fiducial argu- 
ment solves this equation for the unknown parameter fi to get fj, = y — Z; 
for example, suppose y = 4.8, then fj, = 4.8 — Z would suggest \x ~ iV(4.8, 1). 
While the actual value of Z is unknown, the distribution of Z is fully known 
and can be used to frame a distribution on the unknown parameter \x. This 
distribution on \x is known as the fiducial distribution. 

The generalized fiducial recipe starts with a data-generating equation, also 
referred to as the structural equation, which defines the relationship between 
the data and the parameters. Let Y be a random vector indexed by param- 
eter (s) £ € S, and then assume Y can be represented as Y = G(£, U), where 
G is a jointly measurable function indicating the structural equation, and U 
is a random element with a fully known distribution (void of unknown pa- 
rameters). In this paper, the function G will take the form of a normal linear 
mixed model, and the random components U will be standard normal ran- 
dom variables; see equation (2.1). Following the fiducial argument, we define 
a set-valued function, the inverse image of G, as Q(y, u) = {£ : y = G(£, u)}, 
where y is the observed data, and u is an arbitrary realization of U. The 
set-function Q(y, u) is then used to define the fiducial distribution on the pa- 
rameter space. Since the distribution of U is completely known, independent 
copies of U, U* , can be generated to produce a random sample of Q(y, u*) 
for the given data y (where u* is a realization of U*). There are several 
sources of nonuniqueness in this framework. In particular, nonuniqueness 
could occur if Q has more than one element, if Q is empty, or due to the 
definition of the structural equation. Nonuniqueness due to the definition of 
the structural equation will not be addressed here as we assume the form of 
the model is known (i.e., normal linear mixed model). To resolve the case 
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when there is more than one element in Q, we can define a rule, call it V, 
for selecting an element of Q. Furthermore, since the parameters £ are fixed 
but unknown, there must be some realization of the random variable U such 
that y = G(£,u) has occurred [i.e., {Q(y, u) ^ 0}]. The generalized fiducial 
distribution of £ is defined as 

(1-2) V(Q(y,U*))\{Q(y,U*)^0}. 

Obtaining a random sample from the generalized fiducial distribution as 
defined in (1.2) where the structural equation, G, takes the form of a normal 
linear mixed model is the focus of the proposed algorithm. 

Defining the generalized fiducial distribution as (1.2) leads to a potential 
source of nonuniqueness due to conditioning on events with zero probability 
[i.e., if P({Q(y, u) ^ 0}) = 0]. This is known as the Borel Paradox [Casella 
and Berger (2002)]. Fortunately this can be resolved by noting that most 
data has some degree of known uncertainty due, for example, to the reso- 
lution of the instrument collecting the data or computer storage. Because 
of this, instead of considering the value of a datum, an interval around the 
value can be used [Hannig (2012), Hannig, Iyer and Wang (2007)]. For ex- 
ample, suppose the datum value is y = 1.632 meters measuring the height 
of a woman. If the resolution of the instrument used to measure the woman 
is 0.001 m (i.e., 1 mm), then her actual height is between 1.631 meters and 
1.632 meters (or between 1.632 meters and 1.633 meters, depending on the 
practice of the measurer). 

By considering interval data, the issue of nonuniqueness due to the Borel 
Paradox is resolved since the probability of observing our data will never 
be zero since P(Q((a,b),U*) ^ 0) > P(Y E (a,b)) > where a < b are the 
endpoints of the interval. 

Interval data is not explicitly required for generalized fiducial inference, 
but is useful in the proposed setting of normal linear mixed models. General- 
ized fiducial inference has a number of other applications in various settings 
such as wavelet regression [Hannig and Lee (2009)], confidence intervals 
for extremes [Wandler and Hannig (2012)], metrology [Hannig, Iyer and 
Wang (2007)] and variance component models [E, Hannig and Iyer (2008)], 
which applies the generalized fiducial framework to unbalanced normal lin- 
ear mixed models with two variance components. 

1.2. Sequential Monte Carlo. When integrals of interest are very com- 
plex or unsolvable by analytical methods, simulation-based methods can be 
used. SMC, or particle filters, is a collection of simulation methods used to 
sample from an evolving target distribution (i.e., the distribution of interest) 
accomplished by propagating a system of weighted particles through time or 
some other index. A solid introduction of and applications to SMC methods 
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can be found in Doucet, de Preitas and Gordon (2001). There are many di- 
mensions to the theory and applications of SMC algorithms [Chopin (2002, 
2004), Del Moral, Doucet and Jasra (2006), Douc and Moulines (2008), 
Kong, Liu and Wong (1994), Liu and Chen (1998), Liu and West (2001)], 
but a simplified introduction is presented below. 

Suppose one desires to make inferences about some population based on 
data Y. A particle system {z^ J ' ,w^} for J = 1, . . . , N particles is a collec- 
tion of N weighted random variables (with weights w^) such that 

lim 71 ) ^EMz)}= l(z)d7r(z), 

AT— >oo 22 j=l w[ ' J 

where tt is the target distribution, and 7 is some measurable function, when 
this expectation exists. Since it is often difficult to sample directly from 
the target distribution n, it becomes necessary to find some proposal dis- 
tribution tt to sample the particles. The weights w( J > are determined in 
order to re- weight the sampled particles back to the target density [i.e., 
= ir(z^) /tt(z^)]. This is the general idea of importance sampling 
(IS). If the target distribution is evolving with some time index, an iterative 
approach to calculating the weights is performed. This is known as sequen- 
tial importance sampling (SIS). Unfortunately, with each time step, more 
variation is incorporated into the particle system, and the weights degener- 
ate, leaving most particles with weights of essentially zero. The degeneracy 
is often measured by the effective sample size (ESS), which is a measure of 
the distribution of the weights of the particles. Kong, Liu and Wong (1994) 
presents the ESS as having an inverse relation with the coefficient of vari- 
ation of the particle weights, and proved that this coefficient of variation 
increases as the time index increases (i.e., as more data becomes available) 
in the SIS setting. An intuitive explanation of ESS can be found in Liu and 
Chen (1995). 

SMC builds on ideas of sequential importance sampling (SIS) by incorpo- 
rating a resampling step to resolve issues with the degeneracy of the parti- 
cle system. Once the ESS for the particle system has dropped below some 
designated threshold or at some pre-specified time, the particle system is 
resampled, removing inefficient particles with low weights and replicating 
the particles with higher weights [Liu and Chen (1995)]. There are various 
methods for resampling with the most basic being multinomial resampling, 
which resamples particles based on the normalized importance weights; see 
Douc, Cappe and Moulines (2005) for a comparison of several resampling 
methods. 

Examples of general SMC algorithms can be found in Del Moral, Doucet 
and Jasra (2006) or Jasra, Stephens and Holmes (2007). The main idea of 
SMC methods is to iteratively target a sequence of distributions {ftt}t£Z+> 
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where 7Tj is often some distribution based on the data available up to time t. 
The algorithm comprises three main sections after the initialization step: 
sampling, correction and resampling. The sampling step arises at a new 
time step t when particles are sampled from some evolving proposal dis- 
tribution 7Tj. The correction step is concerned with the calculation of the 
weights and the idea of reweighting the particles to target the desired dis- 
tribution at time t, TTt- The resampling step is performed when the ESS of 
the particle system falls below some desired threshold T (e.g., T = N/2). 
The asymptotic correctness for SMC algorithms can be found in Douc and 
Moulines (2008). 

2. Method. 

2.1. Introduction. The form of the normal linear mixed model from equa- 
tion (1.1) is adapted to work in the generalized fiducial inference setting as 

r h 

(2.1) Y = X/3 + ^T^V-,^, 

i=i j=i 

where X is a known nxp fixed-effects design matrix, j3 is the px 1 vector of 
fixed effects, Vij is the n X 1 design vector for level j of random effect i, li is 
the number of levels per random effect i, af is the variance of random effect 
i and the Zij are independent and identically distributed standard normal 
random variables. We will derive a framework for inference on the unknown 
parameters (/3 and <7j, i = 1, . . . , r) of this model, which will be applicable to 
both the balanced and the unbalanced case. (The design is balanced if there 
is an equal number of observations in each level of each effect, otherwise the 
design is unbalanced.) In addition, the covariance matrices for the random 
components (G and R above) are identity matrices with a rank equal to 
the number of levels of its corresponding random effect with l r = n. The 
Vj, i = 1, . .. , r, allow for coefficients for the random effects, and correlation 
structure can be incorporated into the model by including additional effects 
with design matrices that account for the correlation. We note that these 
additional assumptions are related to the implementation of the proposed 
algorithm, and are not restrictions of the generalized fiducial framework. 

Consider the following example illustrating the connection between (1.1) 
and (2.1) in the case of a one-way random effects model. A one-way random 
effects model is conventionally written as 

(2.2) yij = fi + a>i + £ij, i = l,...,a,j = l,...,nj, 

with unknown mean //, random effect a ~ N(0, cr^) where is unknown, a is 
the number of levels of a, rii is the number of observations in level i and error 
terms ~ N(0, a^) where is also unknown, and a and e are independent 
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[Jiang (2007)]. This can be structured in the form of equation (2.1) as 

l\ n 
Y = X/3 + £71 V h3 Z l,3 + °1 V 2J Z 2,j, 
3=1 3=1 

where f3 = fj, is the overall mean, X = l n (an n x 1 vector of ones), l\ = a 
is the number of levels for the first random effect and V\j indicates which 
observations are in level j with random effect variance a\ = a\ . The second 
random effect corresponds to the error, and hence Vi,. = In with c\ as the 
error variance component. The z\ . and Z2 t . are the i.i.d. standard normal 
random variables. 

The SMC algorithm presented in this section is seeking a weighted sample 
of particles {-^ip, W^ J t (where w[ J ^ is the unnormalized importance 

weight for particle z["? t ) from the generalized fiducial distribution of the 
unknown parameters in the normal linear mixed model. Once this sample 
of N weighted particles is obtained, inference procedures such as confidence 
intervals and parameter estimates can be performed on any of the unknown 
parameters or functions of parameters. For example, parameter estimates 
can be determined by taking a weighted average of the particles with the as- 
sociated (normalized) weights. A C% confidence interval can be found easily 
for each parameter by ordering the particles and finding the particle values 
9l and 9jj such that the sum of the normalized weights for the particles 
between 9i and 9{j is C%. 

2.2. Algorithm. The algorithm to obtain a generalized fiducial sample for 
the unknown parameters of normal linear mixed models is outlined below. 
As discussed earlier, the data Y are not observed exactly, but rather intervals 
around the data are determined by the level of uncertainty of the measure- 
ment (e.g. due to the resolution of the instrument used). The structural 
equation formed as interval data for t = 1 , . . . , n with i = 1, . . . , r random 
effects is 

r h 

(2.3) a t <Y t = X t l3 + J2^Yl 

i=l j=l 

where the random effect design vector component Vij t indicates the jih 
level of a random effect i for the tth element of the data vector, Xt is the 
tth row of the fixed effect design matrix X and normal random 

variable for level j of random effect i. Each datum can have one or more 
random components, and so we write Z\-& = (Z\, . . . , Zt) with capital letters 
to indicate the possible vector nature of each Z^ for k = 1, . . . ,t. In the case 
that r > 1, Z\- t is vectorized to be a Yli=i h x 1 vector. Also, for notational 
convenience, Z^ will represent all Zij for i = 1, . . . , r and j = 1, . . . , Zj that 
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are not present or shared with any datum less than k, which will always 
at least include the error effect, denoted z k ^ T . It will be necessary at times 
to reference the nonerror random effects, and they will be denoted Zk i :r —i, 
representing all the nonerror random effects first introduced at time k. The 
goal is to generate a sample of the Zi j for i = 1, . . . , r and j = 1, . . . , Zj such 
that at <Yt <bt for t = 1, . . . , n. 

With a sample of size n, the generalized fiducial distribution on the pa- 
rameter space can be described as 

(2.4) V(Q((a, b] 1:n , Z 1:n )) | {Q((a, h] 1:n ,Z 1:n ) + 0}, 

where we define the set function Q as the set containing the values of the 
parameters that satisfy equation (2.3), given the data and random compo- 
nent Z. Generating a sample from (2.4) is equivalent to generating the Z 
such that the set Q is nonempty, and this results in a target distribution at 
time t written as 

(2.5) 7ri !t (Zi :t | (a, b] 1:4 ) = ir 1:t (Z ltt ) cx exp(-(Zf :t ■ Z 1:t )/2) ■ l Ct {Zv.t), 

where Ic t (") i s an indicator random variable for the set Ct, where Ct is the 
set of Z\-t such that Qt is not empty. This is equivalent to Cf = {Z\- t : 3/3, <7j 
so that a k < X k f3 + J2l=i a i Y?j=i v i,j,k z i,i <h,k = l,...,t}. The restriction 

that is nonempty can be translated into truncating the possible val- 
ues of the particle corresponding to the error random effect to the interval 
defined by 

mj7 (J) 7 (J) s 

That is, mt(z['^_ 1 , Z\ f-"^?, x ) and M t (z[ J ^_ 1 , Z^. r _ 1 ) are the minimum and 

maximum possible values of . 

The proposal distribution used is the standard Cauchy distribution due to 
improved computational stability of sampling in the tails over the more nat- 
ural choice of a standard normal distribution. Specifically, the z[ J J are sam- 
pled from a standard Cauchy truncated to (?n t (z[ J ^_ 1 , zj. J 1 ) . r _ 1 ), M t (z[ J ^_ l , 
^ti-r— i)) f° r t > P + r (i-e., for t greater than the dimension of the parameter 
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space); otherwise z\ T , like Z t -±. r _±, is sampled from a standard normal dis- 
tribution. The conditional proposal density at time t for t > p + r is defined 

as 

^t\l:t-l( Z t I ^l:t-l,(a,b]l;t) 
= ^t\l:t-l{ Z t) 

CXexp(-(2' t T 1:r _ 1 • Z t ,\; T -l)/2) X I( mt (Z 1;t _ 1 ,Z i , 1:r _ 1 ),Af t (Z 1:t _ 1 ,Z f , 1:r „ 1 ))(^,r) 

x [(1 + zf^iFiMtiZt^Zt^t)) - F(m t (Z 1:t -i,Z t , 1:r - 1 )))]- 1 , 

where F is the standard Cauchy cumulative distribution function. Then the 
full proposal density at time t is 

t 

(2.6) TTi:t(Z 1:t | (a,b]i :t ) =7r 1 :t(Zi :t ) = 7ri (Zj | (a, b]i :i ). 

The weights are defined as the ratio of the full joint target density to 
the full joint proposal density at time t. More specifically, the weights are 
derived as 

(2.7) W 1 ..t=Tr 1 .. t /n 1:f . 

The resulting sequential updating factor to Wx-.t-i is Wt = exp(— r /2)(l + 
z^ r )(F(M t (Z 1:t -i, Z tl i:r-i)) - F(mt{Z\-t~i,Z t ,\;r-\))) where F is the stan- 
dard Cauchy cumulative distribution functions. More details regarding the 
derivation of these weights can be found in Appendix C. 

Standard SMC resampling finds particles according to the distribution of 
weights at a given time step, copies the particle and then assigns the resam- 
pled particles equal weight. By copying particles in this setting, we would 
not end up with an appropriate distribution on the parameter space. Intu- 
itively, this is because after each time step, each particle implies the set, or 
geometrically the polyhedron, of possible values of the unknown parameters 
given the generated Z^f t . 

If the particles are simply copied, the distribution of polyhedrons will be 
concentrated in a few areas due to particles with initially higher weight, and 
will not be able to move from those regions because subsequent particles 
would continue to define subsets of the copied polyhedrons, as outlined in 
the algorithm presented above. Hence rather than copy the selected particles 
exactly, we alter them in a certain way in order to retain properties of heavy 
particles, while still allowing for an appropriate sample of {Qt }j=v It 
can be thought of as a Gibbs-sampling step in a noncoordinate direction 
determined by the selected particle. The precise mathematics of this step 
can be found in Appendix A. 
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The proposed algorithm targets the generalized fiducial distribution of 
the unknown parameters of a normal linear mixed model of (2.4) displayed 
in (2.5). The following theorem confirms that the weighted particle system 
from the proposed algorithm achieves this as the sample size N approaches 
infinity. The proof is in Appendix C. 

Theorem 2.1. Given a weighted sample {z[ J ^, }jLi obtained using 
the algorithm presented above targeting (2.5), then for any bounded, mea- 
surable function f as N — > oo, 



This result holds for slightly weaker conditions, which are outlined in 
Appendix C. When the data are i.i.d. (e.g., when the error effect is the 
only random component), the confidence intervals based on the generalized 
fiducial distribution are asymptotically correct [Hannig (2012)]. When the 
data are not i.i.d., previous experience and simulation results suggest that 
the generalized fiducial method presented above still leads to asymptotically 
correct inference as the sample size n increases; see the supplemental docu- 
ment [Cisewski and Hannig (2012)] for a short simulation study investigating 
asymptotic properties of the proposed method and algorithm. The asymp- 
totic exactness of generalized fiducial intervals for two-component normal 
linear mixed models was established in E, Hannig and Iyer (2008); asymp- 
totic exactness of generalized fiducial intervals for normal linear mixed mod- 
els is a topic of future research. 

3. Simulation study and applications. This simulation study has two 
parts. In the first part, we consider the unbalanced two- fold nested model 
with model designs selected from Hernandez, Burdick and Birch (1992). In 
the second part, we use the unbalanced two-factor crossed design with an 
interaction term with designs selected from Hernandez and Burdick (1993); 
both sets of designs include varying levels of imbalance. In addition to the 
classical, ANOVA-based methods proposed in Hernandez, Burdick and Birch 
(1992) and Hernandez and Burdick (1993), we compare the performance of 
our method to the /i-likelihood approach of Lee and Nelder (1996), and a 
Bayesian method proposed in Gelman (2006). The purpose of this study is 
to compare the small-sample performance of the proposed method with cur- 
rent methods using models with varying levels of imbalance. The methods 
were compared using frequentist repeated sampling properties. Specifically, 
performance will be compared based on empirical coverage of confidence 
intervals and average confidence interval length. It is understood that the 
selection of a prior distribution influences the behavior of the posterior; the 
priors were selected based on recommendations in the literature for nor- 
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mal linear mixed models as noted above. While Bayesian methods do not 
necessarily maintain frequentist properties, many practitioners interpret re- 
sults from Bayesian analyses as approximately frequentist (i.e., they expect 
repeated-sampling properties to approximately hold) due to the Bernstein- 
von Mises theorem [Le Cam (1986), van der Vaart (2007)], and so performing 
well in a frequentist sense has appeal. There are a number of examples inves- 
tigating frequentist performance of Bayesian methodology such as Diaconis 
and Freedman (1986a, 1986b), Ghosal, Ghosh and van der Vaart (2000) and 
Mossman and Berger (2001). 

It is important to note that the proposed method is not restricted to the 
model designs selected for this study, and can be applied to any normal linear 
mixed model that satisfies the assumptions from previous sections, while the 
included ANOVA methods were developed specifically for the model designs 
used in this study. A more efficient algorithm than the proposed method may 
be possible for specific model designs, but one of our goals was to present a 
mode of inference for any normal linear mixed model design. 

As presented below, the proposed method tends to be conservative with 
comparable or shorter intervals than the nonfiducial methods used in the 
study. 

3.1. Unbalanced two-fold nested model. For the first part of the simula- 
tion study, we consider the unbalanced two-fold nested linear model 

(3.1) yijk = fJ> + ati + Pij + £ijk 

for i = 1, . . . , I, j = 1, . . . , Ji, and k = 1, . . . , K^, where ft is an unknown 
constant and a\ ~ N(0, a^), /%j ~ A^(0, cr|) and ~ A^(0, a 2 e ). 

Table 1 displays the model designs used in this part of the simulation 
study. Five model designs of Hernandez, Burdick and Birch (1992) were se- 
lected to cover different levels of imbalance both in the number of nested 
groups (Ji) and the number of observations within each group (Kij). The 



Table 1 

Model designs used in the two-fold nested model of (3.1) 



Design 




4>2 




I 


Ji 


Kij 


n 


MI-1 


0.9000 


0.8889 


0.8090 


5 


2,1,1,1,1 


4,4,2,2,2,2 


16 


MI-2 


0.7778 


0.7337 


0.6076 


3 


4,2,1 


1,5,5,5,1,5,1 


23 


MI-3 


1.0000 


1.0000 


1.0000 


3 


3,3,3 


2,2,2,2,2,2,2,2,2 


18 


MI-4 


0.4444 


1.0000 


0.4444 


6 


1,1,1,1,1,7 


2,2,2,2,2,2 


24 














2,2,2,2,2,2 




MI-5 


1.0000 


0.4444 


0.4444 


3 


2,2,2 


1,1,1,1,1,7 


12 



Note: 4>\ and <f>2 reflect the degree of imbalance due to Ji and Kij , respectively, and <j> is 
an overall measure of imbalance. See (3.1) for definitions of /, Ji and Kij; note sample 
size (n) = Y,iY,j K ij- 
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parameters <f>\ and <j>2 reflect the degree of imbalance due to J{ and Kij, 
respectively. The measures of imbalance listed is based on methods pre- 
sented in Khuri (1987) where values range from to 1, and smaller val- 
ues suggest a greater degree of imbalance. The parameters' values used 
in this part of the study are /z = 0, and the following combinations of 
(o£,o$,of): PI-1 = (0.2,0.1,0.7), PI-2 = (0.4, 0.3, 0.3), PI-3 = (0.2, 0.7, 0.1), 
PI-4 = (25, 4, 16) and PI-5 = (1, 1, 1). 

For each model and parameter design combination, 2000 independent 
data sets were generated, and 5000 particles were simulated for the pro- 
posed method. Hernandez, Burdick and Birch (1992) present two methods 
for determining confidence intervals for u^, and three methods for confi- 
dence intervals on cr? . One of the methods is based on the confidence interval 
construction presented in Ting et al. (1990) for balanced designs (denoted 
TYPEI). The other method invokes unweighted sum of squares and is de- 
noted USS. We do not consider the third method presented in Hernandez, 
Burdick and Birch (1992) for confidence intervals on <x| because there is 

not an analogous method for a 2 a . We note that for unbalanced designs, the 
decomposition of the sum-of-squares is not unique and the desired distribu- 
tional properties (independence and chi-squared) do not generally hold. 

The /i-likelihood approach of Lee and Nelder (1996) was implemented us- 
ing the R package hglm, and the results will be referenced as HLMM. The 
/i-likelihood methodology is an approach for extending the likelihood in the 
case of additional, unobserved, random components [Lee, Nelder and Paw- 
itan (2006)]. In the R package hglm, inference on the variance components 
is performed on the log scale. This package was selected because it allows 
for multiple random effects terms, and it includes standard errors on the 
estimates of variance components. 

A Bayesian method is also considered for comparison. Bayesian hierarchi- 
cal models provide a means of constructing confidence intervals for random- 
effects models. Part of the art of the Bayesian methodology is in selecting 
appropriate prior distributions for the unknown parameters. For inference on 
the unknown variance component parameters when there is no prior informa- 
tion available (i.e., when seeking a noninformative prior), Gelman (2006) rec- 
ommends employing a uniform prior distribution on the standard deviation 
parameters when there are a sufficient number of groups (at least 5); other- 
wise, a half-t distribution is suggested. Per the recommendation of Gelman 
(2006), both uniform and half-t priors are considered (denoted BAYI1.5 and 
BAYI3, and BAY2i.5 and BAY23, respectively, where the subscripts 1.5 and 
3 specifies the prior scale variable as explained in Appendix B). The R pack- 
age rjags is used to implement this method; see Appendix B for more details. 

Performance is based on the empirical coverage of (1 — a) 100% confi- 
dence intervals and average interval length for the parameters of inter- 
est 6. We define a lower-tailed (1 — a) 100% confidence interval on 9 as 
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the interval (—00, U a ] such that P(— 00 < 9 < U a ) = 1 — a, an upper-tailed 
(1 — a) 100% confidence interval on 9 as the interval [L a ,oo) such that 
P(L a <9< 00) = 1 — a and a two-sided equal-tailed confidence interval on 
9 as the interval [L a / 2 , U a j 2 ] such that P(L a / 2 <9< U a / 2 ) = I — a. Based 
on the normal approximation to the binomial distribution, we will consider 
empirical coverage between 94% and 96% appropriate for 95% two-sided 
confidence intervals. 

A summary of the two-sided 95% confidence interval results are displayed 
in Figure 1. In addition, the supplemental document [Cisewski and Hannig 
(2012)] includes figures with the results summarized by parameter along with 
the complete raw data results. Average interval lengths are not included in 
the displayed results for HLMM because the excessive lengths would skew 
the scale of the plots; however, the average lengths are displayed in the raw 
data results in the supplemental document [Cisewski and Hannig (2012)]. 

HLMM has empirical coverage below the stated level, and the overall 
longest average interval lengths. BAYI1.5, BAYI3, BAY2i.5 and BAY23 tend 
to be conservative with the longest average interval lengths after HLMM. 
USS and TYPEI maintain the stated coverage well, and FID tends to be 
conservative. Even though FID is conservative, its average interval lengths 
are comparable or shorter than those of USS and TYPEI. However, the 
average interval lengths of USS and TYPEI for cj| are surprisingly wider 
than BAYI1.5, BAYI3, BAY2 L5 , BAY2 3 and FID for MI-1, as revealed in 
the plot of the average lengths in Figure 1(a). This is due to the model 
design; specifically, the derivation of the confidence intervals results in the 
degrees of freedom of 1 for the nested factor (calculated as ^i=i Ji ~ I = !)• 
BAYI1.5, BAYI3, BAY2i.5, BAY2 3 and FID do not appear to have this 
issue. Upper and lower one-sided confidence interval results for a\ and a\ 
for USS and TYPEI tend to stay within the stated level of coverage, and 
BAYI1.5, BAYI3, BAY2 L5 , BAY2 3 and FID range from staying within the 
stated level of coverage to very conservative. 

The proposed method, while maintaining conservative coverage, has av- 
erage interval lengths that are competitive or better than the other methods 
used in this part of the study. The proposed method offers an easily gen- 
eralizable framework and provides intervals for fixed effects and the error 
variance component, unlike the methods presented in Hernandez, Burdick 
and Birch (1992). While the conservative coverage for the Bayesian methods 
can be deemed acceptable, their average interval lengths tend to be wider 
than the proposed method. 

3.1.1. Application 1. In addition to the simulation study, we consider 
the application of model (3.1) presented in Hernandez, Burdick and Birch 
(1992) concerning the blood pH of female mice offspring. Fifteen dams were 
mated with 2 or 3 sires, where each sire was only used for one dam (i.e., 
37 sires were used in the experiment), and the purpose of the study was 



(a) (b) 




USS TYPEI BAY1J1 5} BAY1_3 BAY2J1.5} BAY2_3 HB BAY1J1.5} BAY13 BAY2J1.5} BAY23 

Method Method 

Fig. 1. Summary of simulation results, (a) Displays the combined results for a\ and cr| for the two-fold nested model, (b) Displays the 
combined results for o\, cr^ and a 2 a p for the two-factor crossed design with interaction. The first row is the empirical coverage probabilities 
for 95% two-sided confidence intervals. The second row, first column is the log 10 of the average interval lengths divided by the average 
interval lengths of FID, and the second row, second column is the average interval lengths divided by the average interval lengths of FID. 
Each value used in a box plot corresponds to a particular model design and parameter combination for each nonerror variance component. 
The average interval lengths for HLMM were not included in the displayed results due to their excessive lengths. All the simulation results 
are available in the supplemental document [Cisewski and Hannig (2012)]. 
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Table 2 

Two-fold nested model: real data example 



V (11 • CUII1JJ. 


IVltitllULl 




i-OlUtJU / dvc. It311> 


UJJjJtJl /lUWcl 


Co 


uss 


(2.30, 28.56) 


0.953/25.1 


0.949/0.958 




TYPEI 


(1.94, 26.23) 


0.950/25.2 


0.949/0.957 




BAYli.s 


(1.73, 30.72) 


0.955/29.2 


0.948/0.959 




BAY1 3 


(1.51, 30.21) 


0.956/37.2 


0.948/0.962 




BAY2i.5 


(1.56, 30.02) 


0.955/27.3 


0.948/0.960 




BAY2 3 


(1.76, 30.04) 


0.955/27.5 


0.950/0.961 




FID 


(1.53, 26.67) 


0.947/24.5 


0.958/0.947 




USS 


(0.00, 11.56) 


0.961/10.9 


0.952/0.952 




TYPEI 


(0.00, 11.26) 


0.964/12.4 


0.953/0.953 




BAYI1.5 


(0.17, 11.81) 


0.976/11.2 


0.956/0.994 




BAYI3 


(0.04, 12.55) 


0.980/11.3 


0.959/0.996 




BAY2i.s 


(0.01, 12.49) 


0.982/11.3 


0.958/0.994 




BAY2 3 


(0.01, 11.92) 


0.983/11.3 


0.958/0.995 




FID 


(0.19, 10.54) 


0.974/10.7 


0.951/0.986 



Note: The 95% intervals are based on the actual data while the remaining information 
are the empirical results from 2000 independently generated data set using the REML 
estimates for each parameter. The results are the empirical coverage and average interval 
lengths of 95% confidence intervals. 



to determine if the variability in the blood pH of the female offspring is, in 
part, due to the variability in the mother. There is imbalance in the data due 
to the number of sires mated with each dam (2 or 3); also note the natural 
imbalance in the data resulting from the number of female offspring. 

The 95% confidence intervals based on the real data are presented in Ta- 
ble 2. An example of the generalized fiducial distribution for a\ is displayed 
in Figure 2. This highlights one of the advantages of the proposed method 
over classical methods (and shared with Bayesian methods), which is a dis- 
tribution on the parameter space allowing for inferences similar to those 
made using Bayesian posterior distributions. 

In order to evaluate the empirical coverage of the proposed method, we 
perform a simulation study using the REML estimates for all the param- 
eters (/j = 44.92, a 2 a = 8.90, aj = 2.65 and of = 24.81). Simulating 2000 
independent data sets with the noted parameter values, we find the empir- 
ical coverage using USS, TYPEI, BAY1 L5 , BAYI3, BAY2i. 5 , BAY2 3 and 
FID, and the average lengths of the two-sided intervals. The results of the 
simulation study are also found in Table 2. 

The confidence interval coverage and average lengths for a\ are compara- 
ble for all methods with BAYI3 having the longest average interval length. 
For oj, BAYI1.5, BAYI3, BAY2 L5 , BAY2 3 and FID tend to be more conser- 
vative while USS and TYPEI correctly maintain the stated level of coverage. 
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2 

Fig. 2. The generalized fiducial distribution of 5000 generated particles for a\ using the 
real data described in Section 3.1.1. The solid line is the normal kernel density estimate 
of the distribution with a point mass at zero. The two vertical dashed lines represent the 
lower and upper bounds for 95% confidence intervals based on the weights of the generated 
particles. 

USS, TYPEI and FID have comparable interval lengths, and the average 
lengths of BAYI1.5, BAY1 3 , BAY2 L5 , BAY2 3 are slightly longer. 

3.2. Unbalanced two-factor crossed design with interaction. In this part 
of the simulation study, we consider the unbalanced two-factor crossed de- 
signs with interaction written as 

(3.2) Yij k = /i + «i + Pj + {ap)ij + e ijk 

for i = 1, . . . , I, j = 1, . . . , J and k = 1, . . . , Kij, where fi is an unknown 
constant and a.{ ~ N(0,a^), ~ 2V(0,ct|), (a/3)ij ~ iV(0, cr^) and e^k ~ 
N(0,o-*). 

This model is presented in Hernandez and Burdick (1993), where the 
authors propose a method based on unweighted sum of squares to construct 
confidence intervals for cr^, cr| and a^a- The method they propose is based 

on intervals for balanced designs presented by Ting et al. (1990). In the 
simulation study, this method will be called HB. As with (3.1), HLMM 
from Lee and Nelder (1996), and BAY1 L5 , BAY1 3 , BAY2i. 5 and BAY2 3 
from Gelman (2006) will be used as a comparison. 

Table 3 displays the model designs used in this part of the study, and, 
again, the overall measure of imbalance (0) proposed in Khuri (1987) is 
displayed for each design. The parameters values used in this part of the 
study are fi = 0, and the following combinations of (<7„, aj, cr^, a^): PII-1 = 
(0.1,0.5,0.1,0.3), PII-2 = (0.1,0.3,0.1,0.5), PII-3 = (0.1, 0.1, 0.3, 0.5), PII-4 = 
(0.1,0.1,0.5,0.3), and PII-5 = (1, 1, 1, 1). 
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Table 3 

Model designs used in the two- factor crossed design with interaction model of (3.2) 



Design 




I 


J 




n 


MII-1 


0.8768 


4 


3 


2,1,3/2,1,1/2,2,2/1,2,3 


22 


MII-2 


0.6667 


3 


3 


4,1,1/4,1,1/4,1,1 


18 


MII-3 


0.6667 


3 


3 


4,4,4/1,1,1/1,1,1 


18 


MII-4 


0.4011 


3 


4 


8,1,1,1/1,1,1,1/1,1,1,1 


19 


MII-5 


0.7619 


5 


3 


1,2,2/5,2,7/2,2,2/2,4,2/3,2,2 


40 


MII-6 


1.000 


3 


3 


2,2,2/2,2,2/2,2,2 


18 



Note: The parameter <j> is an overall measure of imbalance of the model. See (3.2) for 
definitions of /, J and K^; note sample size (n) = YJ 4 YJ^ Kij. 



For each design and set of parameter values, 2000 independent data sets 
were generated, and 5000 particles were simulated for the proposed method. 
As before, performance is based on the empirical coverage of (1 — a) 100% 
confidence intervals and average interval length for the parameter of inter- 
est 9. Based on the normal approximation to the binomial distribution, we 
will consider empirical coverage between 94% and 96% appropriate for 95% 
two-sided confidence intervals. 

A summary of the two-sided 95% confidence interval results for a^, cr| 
and o 2 a p are displayed in Figure 1. In addition, the supplemental document 
[Cisewski and Hannig (2012)] includes figures with the results summarized 
by parameter along with the complete raw data results. Average interval 
lengths are not included in the displayed results for HLMM because the 
excessive lengths would skew the scale of the plots; however, the average 
lengths are displayed in the raw data results in the supplemental document 
[Cisewski and Hannig (2012)]. 

HLMM has empirical coverage below the stated level, and the longest 
average interval lengths. While HB maintains the stated coverage and FID 
tends to be more conservative, they have comparable average interval lengths. 
BAYI1.5, BAYI3, BAY2i,5, and BAY23 are conservative with the longest av- 
erage interval lengths after HLMM for a\ and cr|; while still conservative 
for c^a, the average interval lengths are shorter than HB, but longer than 
FID. One-sided confidence interval results for HB maintains the stated cov- 
erage, and BAYI1.5, BAYI3, BAY2 L5 , BAY2 3 and FID tends to be within 
the stated coverage to very conservative. 

4. Conclusion. Even with the long history of inference procedures for 
normal linear mixed models, a good-performing, unified inference method is 
lacking. ANOVA-based methods offer, what tend to be, model-specific so- 
lutions. While Bayesian methods allow for solutions to very complex mod- 
els, determining an appropriate prior distribution can be confusing for the 
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nonstatistician practitioner. In addition, for the models considered in the 
simulation study and the prior selected based on recommendations in the 
literature, the Bayesian interval lengths were not generally competitive with 
the other methods used in the study. The proposed method allows for con- 
fidence interval estimation for all parameters of balanced and unbalanced 
normal linear mixed models. In general, our recommendation is to use the 
proposed method because of its apparent robustness to design imbalance, 
good small sample properties and flexibility of inference due to a fiducial 
distribution on the parameter space. If the design is balanced and only 
confidence intervals are desired, an ANOVA-based method would provide a 
computationally efficient solution. 

It is interesting to note that even though more variation was incorporated 
into the data for the proposed method due to its acknowledgment of known 
uncertainty using intervals, in the simulation study, the proposed method 
tended to have conservative coverage, but the average interval lengths were 
comparable or shorter than the other methods that assumed the data are 
observed exactly. The currently implemented algorithm is suitable for 9 or 
fewer total parameters, but the method does not limit the number of pa- 
rameters. A MATLAB implementation of the proposed algorithm is avail- 
able on the author's website at http://www.unc.edu/~hannig/download/ 
LinearMixedModel_MATLAB . zip. 

APPENDIX A: RESAMPLING ALTERATION STEP 

The particle to be resampled is decomposed into an orthogonal projection 
onto a certain space (discussed below) and the part orthogonal to that space, 
and then the distributional properties that arise from the decomposition are 
used to alter the resampled particle. The alteration step of the proposed algo- 
rithm is performed in such a way that it still solves the system of inequalities 
of (2.3) up to time t using the following idea (to ease the notational complex- 
ity, we do not include the dependence of the variables on t). Suppose particle 
L is selected to be resampled (for an L between 1 and N). For each random 

effect, e, let Y = X'/3' + (rVZ( L \ where X' = [X, {^ l j =1 V^z^}^}, f3' = 

W, Wi}i^e)', (y = o e and vz( L ) = Y!j=i V ed Z ir In order to alter z(i) > we 
first find the basis vectors, rj, for the null space M = null[— X', V] = rj where 
for matrix A = [— X', V], null(A) is the set {rj-.Arj = 0}, and rj = (7/1, 772)^ 
such that (i) X' • r]i = 0, (ii) V • 772 = 0, (iii) r/J • 772 = I (i.e., 772 is orthonormal) 
and (iv) rank(r/) =rank(r]2)- We perform the following decomposition: 

(A.i) z^=nz^ + \\z^-uz^\ 



z(L)-nz( L )\ 

jace M (i.e., I 

Z( L )),and || • || is the L 2 norm. Define C = r% ■ Z^ (so that, ?? 2 ■ C = UZ^), 



where TLZ( L ' is the projection onto the null space M (i.e., nzW = 77-77 
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D= ||Z( L ) -nZ( L )||, andr= Then if ZW is standard normal, 

C ~ A^ e (0,I) where l e is the number of levels of random effect e of Z^ L > to 

be resampled at time t, and D ~ yj Xt-d wnere d = rank (A/"), and C and 

D are independent by design. The alteration of is accomplished by 
sampling new values of C and D (denoted C and D, resp.) according to 
their distributions determined by the decomposition above, and the altered 
particle is 

(A.2) Z = m ■ C + D • r. 

Notice that if Z^ L ' is a standard normal conditioned on Ct, then so is Z, 
and hence the alteration proposed still targets the correct distribution and 
is a Markovian step. 

Furthermore, the set Q { t L) = {(/3' , ct) : a* < X'p'+aVZM < b h i = 1, . . . ,t} 
can be adjusted noting that if (/3',a) solves a, < X'/3' + ctV(7/2 • C + D-r) < b{ 
for i = 1, . . . , t, then (/3, a) can be found such that ctj < X'/3 + crV(ry2 • C + D • 
t) < 6i for i = 1, . . . , t by considering X'/3 + ctVZ = X,5 + crV(?72 • C + rD) = 
X'/3 + a\{j]2 ■ C + rD) = X'/3 + crVZ. Examining the orthogonal parts first, 
ctV(tD) = ctV(tD) implies Vr(crD - 3?D) = implies ct = ct(D/D). The 
relation between /3 and /3' follows from the remaining portion 

X'/3' + ctV(t/ 2 • C) = X'/3 + 5V(?72 • C) implies 

(A-3) 

X'(/3 - f3') + aVr] 2 {C ■ D/D - C) = 0. 
Noting by definition —IL'rji + Vr/2 = 0, then 
(A.4) - XVr?i(C • D/D - C) + ctVt^C • D/D - C) = 0. 

By combining (A.3) and (A.4), we see that /3 - j3' = -ar]i • (C • D/D - C), 
and hence j3 = /3' — ar/i ■ (C - D/D — C). Hence the sets q\ J \z\\ J ^) are easily 
updated to Q\ J \z^). This procedure is repeated for each random effect. 

APPENDIX B: PRIOR DISTRIBUTION SELECTION 

The R package "rjags" was used to implement the Bayesian methods used 
in the simulation study and applications. Gelman (2006) suggests using a 
uniform prior [i.e., U (0, a)] on the standard deviation parameters when there 
are at least 5 groups and explains that fewer than 3 groups results in an 
improper posterior distribution. Calibration is necessary in selecting the 
parameter a in the prior distribution; we use 1.5 and 3 times the range of 
the data (per the recommendation in Gelman [(2006), page 528] to use a 
value that is "high but not off the scale"), which appears reasonable when 
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reviewing the resulting posterior distributions. In the simulation study, the 
results at these scale are denoted by BAYI1.5 and BAYI3, respectively. For 
example, the hierarchical model for the two-fold nested model of (3.1) is 

Yijk ~ N(fi + ai + Pij,o 2 £ ), i = l,...,I,j = l,...,Ji,k = l,..., Kij, 

ai~N(0,a a ), i = l,...,I, /% ~ N(0,ap), j = l,...,Jj, 

<7a~^ r (0, a), ap~U(0,a). 

For the second Bayesian method, a similar hierarchical model is used. In- 
stead of a uniform distribution on the nonerror variance components, a half- 
Cauchy distribution with scale parameter a set as 1.5 or 3 times the range 
of the data (denoted BAY2i.s and BAY23, resp., in the simulation study) is 
used. 



APPENDIX C: PROOF OF THEOREM 

The proof of the convergence of the proposed SMC algorithm follows 
from ideas presented in Douc and Moulines (2008). Theorem 2.1 will fol- 
low from proving the convergence of the generated particles after each stage 
of the algorithm: sampling, resampling and alteration. The development of 
the particle system using the proposed algorithm does not follow the tradi- 
tional SMC algorithm as presented in Douc and Moulines (2008), Section 2. 
A distinction is seen in the formulation of the proposed weights introduced 
in (2.7) and discussed below. 

Using the derivation of the proposal distribution in and above (2.6), 
the target distribution of (2.5) and noting that Ic t (^l:t) = Ict-i(^l:t-l) ' 
I*(Z M:r _i) • Ifa(Zi*-i,Zt,i„-i)Mt(Zi*-i,Zt,i*--i))( z t,r) I where the * in 
I*(^t,i:r— 1) indicates the lack of restriction to a specific set of values for 
Zti;r-i]) the marginal target distribution at time t is 



oc7ri : t_i J exp(-(Zf ■ Z t )/2) 



X I*(-Zt,l:r-l)I(77l t (^l : t-l,^t,l : r-l),Aft(Xi : t-l,^t,l:r-l))( 2 i,'-) <iZ i 
OC 7Ti:t_i • (*(M t (Zi !t _i,Z t|l!r _i)) - *(mt(Z l!t _i,Zt,l:r-l))), 

where F and $ are the standard Cauchy and standard normal cumulative 
distribution functions, respectively, and @i-t is the normalization factor at 
time t. It then follows that the conditional target distribution at time t is 

= 7ri :t (Zl :t )/7fl;t(Zl : t) 
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OCexp(-(Z t T • Z t )/2) ■ I*(Z tj i :r _i) • I(m t (^i : t-i,^t,i : r-i) ) M t (0i:t-i,Zt, l!r -i))(^,: 
/($(M t (Z 1:t _ 1 ,Z til:r _ 1 ))-$(m t (Z 1:t _ 1 ,Z tjl:r _ 1 ))). 
Finally, following the notation just below (2.7), the derivation of the weights 

at time t is w 1:t = &j = ~ « e XP (-^ r / 2 )(i + 4 r)(jF(M<(Zl:t _ 1? 

^,i :r _i)) - F(mt(Zi :t _i, ^t,i:r-i))) • oc Wi • Wi:t_i. This proof will use 
the above formulation and the following notation and definition. 

A particle system is defined as {z[ J ^ , W r 1 ^}j r =1 with z[ J ^ sampled from 
the proposal distribution TT\-t as defined in (2.6) targeting probability mea- 
sure 7ri : t on (@i :t ,J3(@i:t)), and un-normalized weights w[ J ^ as defined 
in (2.7). Let T^m-.t—i be the marginal proposal density at time t as defined 
above equation (2.6), which follows a Cauchy distribution truncated to the 
region R t = (m, t (Z 1:t -i, Z ti i :r ._i), M t (Z 1:t -i, Z t ,i; r -i)) defined by previously 
sampled particles. For notational convenience, let Q t = X)?=i Wff . Define 

two sigma-fields Fq = 0"({^£?}j =1 , ( a , &]i:t) an d ^j-^oV °~({z[*t }i<k<j, 
(a, b]i : t), for J = 1, . . . , iV. Finally, we define proper set .Bt [Douc and Moulines 
(2008), Section 2.1] where B t = {f e L x {Q l:t ^ 1:t ),F{; |/|) G B t ^} and 
F(Zi :t _i,/) = //(Zi :t _i,Zt)I*(^ ) i :r _i)Ifl t (^ )r )($(M t (Zi !t _i,Z tf i !P _i)) - 

^> (m t , Z t; l :r _l) ))7T t | 1:t _! (cZ^t) . 

Definition C.l. Following Definition 1 of Douc and Moulines (2008), a 
weig hted sample {Z$,W$}j =i,...,N is consistent for the probability mea- 
sure iri : t and the proper set B t if, for any / G B t , Q^ 1 X)j=i ^vt fi^vt) ~^ 
I f(Zi lt )T 1 *(dZ 1 4)±ni l t(f), and max^ =1 A 0. 

Two additional conditions on the set B t will be required to guarantee con- 
sistency for the particle system after the alteration step. Let / G Bt, and II 
be the projection matrix onto the null space defined in the description of the 
resampling step of the algorithm found in Appendix A. E[f(z[ J f) | Tj-i] = 

ff(V2 ■ C + D ^}-^b ^c,D =If(V2 ■ C + Dr$)dn d5 4 hf (Z$), 

_ J Z l:i -UZ l:t II 

where C and D are as defined above (A. 2). For hf to be in Bt, f must be 
selected so that the following two conditions hold for any direction Ti- t : 



f{m-C ' + DTi:t) dirx £ 



d/Ki;t < oo, 



(C.l) l\h f (Z 1:t )\d7r 1:t -- 

F(Z 1:t _ 1 ,h f ) = J (J f(r, 2 -C + DT[ :t )dir d>3 \ -l Ct {Zi:t) 



(C.2) 



($(M t ) - ${m t ))n\i:t-i{dZ t ) < oo, 
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where r[. t = fef^fpSB^T^ • Let B t be the set of / G B t such that (C.l) 

and (C.2) hold. Then, Bt C Bt, and we replace Bt with Bt in the definition 
of Bt+i. Finally, since all bounded functions satisfy (C.l) and (C.2), Bt is 
nonempty. 

The goal is to show the particle system generated from the presented 
algorithm is consistent. This requires the particle system after sampling and 
reweighting to be consistent, after resampling to be consistent and after the 
alteration to be consistent. 

After noting (i) for / € B u E(W$ f(Z$) | Tj-i) = W^E^ f(Z$ ) 

Tj-x) = W[^_ x ■ F(Zut-iJ) for J = 1,...,N, with F is defined above 

and Jy-i = cr({^iit-i}jLi, (o,&]l:t) v ^{{Z^}i<k<J-u (a> b]x :t ), and 
(ii) Wt7r t \i:t-i = (^(Mf) — <&(mt))7rt\i : t-i, consistency after sampling and 
reweighting closely follows the proof of Theorem 1 in Douc and Moulines 
(2008). Consistency after resampling follows directly from Theorem 3 in 
Douc and Moulines (2008) and consistency after the alteration step is ad- 
dressed below in Lemma C.l. 

Lemma C.l (Alteration). Assuming the uniformly weighted sample 
{z[ J ^,l}j =l is consistent for {-Ki-t,Bt), then the altered uniformly weighted 
sample {z[ J ^ , ljjL} is consistent for {-Ki- t ,B- 



Proof. Note that the {Z^\l}j f =1 is the altered particle system, while 

1 are the resampled particles. We note that hf is a function of 
Z\ : t because t\-x is a function of Z± : t, and, at times, it will be necessary to 
write = Ti : t(z[ J J). Recall that C and D are defined by a decomposition 
of the original particle selected to be resampled and are independent by 
design. The C and D are the random variables to be resampled according 
to the target distributions of C and D with the r}^ considered fixed so that 

Zv.t = m - C + DT 1:t . 

The lemma will follow once we show 

N N 



AT 1 £ E[f(Z$) I 7>_x] = N- 1 ]T h f (Z$) 
J=i J=i 

— > J hf(Z 1:t )n 1: t(dZ 1:t ) = J f(Zl :t )ir 1:t (dZ 1:t ). 

This is because trivially E[f(z[ J ^) j Tj-i] = f(z[ J t > ), so all that is needed is 
for: 

(i) N^ 1 Ylj=i hf(Z$) — > j hf(Zi; t )iTi;t{dZi;t) and 

(ii) fh f (Z 1:t )diri :t (Z 1:t ) = f f(Z 1:t )7T 1:t (dZ 1:t ). 
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Point (i) follows from (C.l), (C.2), and because f€Bt so that hf&Bt. Now 
we only need to show point (ii) that Jhf(Zi : t)iri : t{dZi : t) = J f (Zi:t)'R'i:t(dZi:t) ' 
This holds because 

J hf(Z 1:t )dni :t = J h*f(n : t) diT T 

J f(m ■ C + Dn-t) dxc,Dj dir T 
= J f(V2 ■ C + Dn-t) d-K C ,D x dn T 

= J f(Z 1:t )dir 1:t = E ni[t [f(Z 1:t )], 

where h*f(r\^) = hf(Zi-t), the equality from line one to line three follows be- 
cause r(Zi-t) = T(Z\-t) and the equality in the third and fourth lines follows 
by Fubini's theorem because C and D are independent of r. □ 

Acknowledgment. The authors would like to thank the reviewers and 
Professor Hari Iyer for their many helpful comments. 

SUPPLEMENTARY MATERIAL 

Additional simulation results (DOI: 10.1214/12-AOS1030SUPP; .pdf). 
The asymptotic stability of the algorithm, with respect to the sample size 
and the particle sample size, was tested, and the simulation results are in- 
cluded in this document. The raw results for the simulation study in Sec- 
tion 3 are also displayed, along with additional summary figures. 
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